Nasz źródłowy dataset zawiera informacje na temat połowu śledzi w Europie, w ciągu ostatnich 60 lat. Zbiór danych zawiera informacje dotyczące długości złowionych śledzi, dostępności planktonu będącego jego pożywieniem oraz informacje na temat samych połowów i rzeczy takich jak temperatura przy powierzchni wody i jej poziom zasolenia w momencie połowu.
Zbiór ten jest podstawą do analizy i określenia głównych przyczyn stopniowego karłowacenia (zmniejszania się długości) śledzi w Europie. Zbiór składa się z następujących atrybutów:
Po analizie danych, na podstawie diagramu korelacji oraz utworzonego regresora, jednoznacznie stwierdzono, że największą przyczyną karłowacenia śledzi jest wzrost temperatury przy powierzchnii wody. Wpływ również miała oscylacja północno-atlantycka, jednak jest bardzo wysokie przypuszczenie że jest to spowodowane tym iż bezpośrednio wpływa ona na temperaturę przy powierzchnii wody.
Nasze tezy znajdują potwierdzenie w wielu artykułach co sugeruje poprawność naszej analizy i modelu.
## X length cfin1 cfin2
## Min. : 0 Min. :19.0 Length:52582 Length:52582
## 1st Qu.:13145 1st Qu.:24.0 Class :character Class :character
## Median :26291 Median :25.5 Mode :character Mode :character
## Mean :26291 Mean :25.3
## 3rd Qu.:39436 3rd Qu.:26.5
## Max. :52581 Max. :32.5
## chel1 chel2 lcop1 lcop2
## Length:52582 Length:52582 Length:52582 Length:52582
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
##
##
##
## fbar recr cumf totaln
## Min. :0.0680 Min. : 140515 Min. :0.06833 Min. : 144137
## 1st Qu.:0.2270 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068
## Median :0.3320 Median : 421391 Median :0.23191 Median : 539558
## Mean :0.3304 Mean : 520367 Mean :0.22981 Mean : 514973
## 3rd Qu.:0.4560 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351
## Max. :0.8490 Max. :1565890 Max. :0.39801 Max. :1015595
## sst sal xmonth nao
## Length:52582 Min. :35.40 Min. : 1.000 Min. :-4.89000
## Class :character 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
## Mode :character Median :35.51 Median : 8.000 Median : 0.20000
## Mean :35.51 Mean : 7.258 Mean :-0.09236
## 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
## Max. :35.61 Max. :12.000 Max. : 5.08000
## 'data.frame': 52582 obs. of 16 variables:
## $ X : int 0 1 2 3 4 5 6 7 8 9 ...
## $ length: num 23 22.5 25 25.5 24 22 24 23.5 22.5 22.5 ...
## $ cfin1 : chr "0.02778" "0.02778" "0.02778" "0.02778" ...
## $ cfin2 : chr "0.27785" "0.27785" "0.27785" "0.27785" ...
## $ chel1 : chr "2.46875" "2.46875" "2.46875" "2.46875" ...
## $ chel2 : chr "?" "21.43548" "21.43548" "21.43548" ...
## $ lcop1 : chr "2.54787" "2.54787" "2.54787" "2.54787" ...
## $ lcop2 : chr "26.35881" "26.35881" "26.35881" "26.35881" ...
## $ fbar : num 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 ...
## $ recr : int 482831 482831 482831 482831 482831 482831 482831 482831 482831 482831 ...
## $ cumf : num 0.306 0.306 0.306 0.306 0.306 ...
## $ totaln: num 267381 267381 267381 267381 267381 ...
## $ sst : chr "14.3069330186" "14.3069330186" "14.3069330186" "14.3069330186" ...
## $ sal : num 35.5 35.5 35.5 35.5 35.5 ...
## $ xmonth: int 7 7 7 7 7 7 7 7 7 7 ...
## $ nao : num 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 ...
Komentarz: Wstępne oględziny zbioru danych sugerują,
że posiada on wartości puste oznaczone symbolem “?” mające symbolizować
błąd/niedostępność/brak danych (NA). Dodatkowo widać, że
niektóre kolumny zawierające liczby są złego typu (character) i powinny
zostać zrzutowane na prawidłowy typ numeryczny.
## X length cfin1 cfin2
## Min. : 0 Min. :19.0 Min. : 0.0000 Min. : 0.0000
## 1st Qu.:13145 1st Qu.:24.0 1st Qu.: 0.0000 1st Qu.: 0.2778
## Median :26291 Median :25.5 Median : 0.1111 Median : 0.7012
## Mean :26291 Mean :25.3 Mean : 0.4458 Mean : 2.0248
## 3rd Qu.:39436 3rd Qu.:26.5 3rd Qu.: 0.3333 3rd Qu.: 1.7936
## Max. :52581 Max. :32.5 Max. :37.6667 Max. :19.3958
## NA's :1581 NA's :1536
## chel1 chel2 lcop1 lcop2
## Min. : 0.000 Min. : 5.238 Min. : 0.3074 Min. : 7.849
## 1st Qu.: 2.469 1st Qu.:13.427 1st Qu.: 2.5479 1st Qu.:17.808
## Median : 5.750 Median :21.673 Median : 7.0000 Median :24.859
## Mean :10.006 Mean :21.221 Mean : 12.8108 Mean :28.419
## 3rd Qu.:11.500 3rd Qu.:27.193 3rd Qu.: 21.2315 3rd Qu.:37.232
## Max. :75.000 Max. :57.706 Max. :115.5833 Max. :68.736
## NA's :1555 NA's :1556 NA's :1653 NA's :1591
## fbar recr cumf totaln
## Min. :0.0680 Min. : 140515 Min. :0.06833 Min. : 144137
## 1st Qu.:0.2270 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068
## Median :0.3320 Median : 421391 Median :0.23191 Median : 539558
## Mean :0.3304 Mean : 520367 Mean :0.22981 Mean : 514973
## 3rd Qu.:0.4560 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351
## Max. :0.8490 Max. :1565890 Max. :0.39801 Max. :1015595
##
## sst sal xmonth nao
## Min. :12.77 Min. :35.40 Min. : 1.000 Min. :-4.89000
## 1st Qu.:13.60 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
## Median :13.86 Median :35.51 Median : 8.000 Median : 0.20000
## Mean :13.87 Mean :35.51 Mean : 7.258 Mean :-0.09236
## 3rd Qu.:14.16 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
## Max. :14.73 Max. :35.61 Max. :12.000 Max. : 5.08000
## NA's :1584
Komentarz: Powyższy fragment kodu zamienia znak “?”
na NA oraz rzutuje na typ numeryczny wartości z kolumn,
które posiadały błędny typ (character). Tak wstępnie wyczyszczone dane
są następnie zapisywane jako nowy dataset o nazwie
tidy_dataset.
Rozkład wartości pustych (NA) w naszym zbiorze
prezentuje się następująco:
## Liczba NA dla danego atrybutu
## X 0
## length 0
## cfin1 1581
## cfin2 1536
## chel1 1555
## chel2 1556
## lcop1 1653
## lcop2 1591
## fbar 0
## recr 0
## cumf 0
## totaln 0
## sst 1584
## sal 0
## xmonth 0
## nao 0
Komentarz: Widać że kolumny zawierające informacje o
dostępności danego planktonu (cfin1,cfin2,chel1,chel2,lcop2,lcop3) oraz
temperatury przy powierzchni wody (sst) zawierają puste wartości. Każdy
atrybut posiadający wartości puste musi zostać przeanalizowany względem
całokształtu zestawu danych by móc określić jakie działanie będzie
najlepsze pod względem statystycznym. Musimy podjąć decyzję czy wartość
NA zostanie np. zastąpiona średnią, ostatnią znaną próbką,
wyzerowana, czy może najlepszym rozwiązaniem będzie usunięcie próbki tj.
powiązanego z nią rekordu.
I tak kolejno dla atrybutów posiadających NA
postanowione zostało:
cfin1,cfin2,chel1,chel2,lcop2,lcop3: oznacza dostępność planktonu określonego gatunku w momencie połowu. Wartość pusta nie może zostać domyślnie zastąpiona zerem ponieważ powstanie sytuacja typu “w jeden dzień zniknął cały plankton, by następnego znowu się pojawić”. Wartość średnia całej kolumny również jest złym pomysłem ponieważ może zawyżyć/obniżyć wartości zebrane w danym oknie czasowym np. wartość średnia atrybutu wyniosła 5, ale w danym oknie czasowym wartość planktonu utrzymywała się między 3 a 4, dlatego wybranie średniej może dać nam mocno przekłamane wyniki.
Najlepszym rozwiązaniem będzie zastąpienie wartości NA
ostatnią wartością niepustą ale z uporządkowanego ciągu
tych wartości, dla każdego atrybutu z osobna. W tym
przypadku uzyskamy sytuację w której najgorszym możliwym przypadkiem
granicznym będzie ostatnia wartość niepusta z poprzedniego analizowanego
połowu z maksymalnym błędem w postaci różnicy dostępności planktonu
między dwoma połowami.
Na przykład dla cfin1 = {15,NA,10}, przy założeniu że
dostępność planktonu liniowo malała, wartość średnia (12.5)
byłaby dobrym rozwiązaniem, natomiast przy takich samych założeniach dla
przypadku:
cfin1 = {15,15,15,15,NA,...,14,14,...,13,...,NA,...,12,...,11,...,NA,...,10},
wybranie wartości średniej dla tego okna czasowego i zestawu próbek
byłoby złym pomysłem ponieważ pierwsze NA powinno mieć
wartość z przedziału <15,14>, a trzecie z przedziału
<11,10>, jedynie dla drugiego NA, które
znajduje się mniej więcej po środku zbioru, wartość średnia (12.5)
spełnia założenie że drugie NA powinno być z przedziału
<13,12>.
sst: oznacza temperaturę przy powierzchnii wody. Ze
względu na mały rozkład wartości <12.77,14.73>
wybranie wartości średniej wydaje się być dobrym pomysłem. Jednakże dla
tego atrybutu zdecydowanie należy wziąć pod uwagę precyzję i czułość.
Zmiana wartości temperatury nawet o 1.0, która mierzona była z dużą
dokładnością może mieć znaczący wpływ na naszą analizę.
Bezpieczniej i lepiej z statystycznego punktu widzenia będzie zastosowanie tego samego podejścia co dla kolumn zawierających informacje o dostępności planktonu. Tutaj w najgorszym możliwym przypadku uzyskamy błąd w postaci różnicy temperatur między dwoma połowami, który powinien być zasadniczo mały jeśli przyjmiemy założenie, że jeśli jeden pomiar był wykonany zimą, to kolejny nie był wykonany dopiero podczas gorącego lata, co w przypadku naszego uporządkowanego zbioru danych jest definitywnie spełnione.
## X length cfin1 cfin2
## Min. : 0 Min. :19.0 Min. : 0.0000 Min. : 0.0000
## 1st Qu.:13145 1st Qu.:24.0 1st Qu.: 0.0000 1st Qu.: 0.2778
## Median :26291 Median :25.5 Median : 0.1333 Median : 0.7012
## Mean :26291 Mean :25.3 Mean : 1.5650 Mean : 2.5323
## 3rd Qu.:39436 3rd Qu.:26.5 3rd Qu.: 0.3603 3rd Qu.: 1.9973
## Max. :52581 Max. :32.5 Max. :37.6667 Max. :19.3958
## chel1 chel2 lcop1 lcop2
## Min. : 0.000 Min. : 5.238 Min. : 0.3074 Min. : 7.849
## 1st Qu.: 2.469 1st Qu.:13.589 1st Qu.: 2.5479 1st Qu.:17.808
## Median : 6.083 Median :21.673 Median : 7.1229 Median :25.338
## Mean :11.928 Mean :22.301 Mean : 16.0416 Mean :29.639
## 3rd Qu.:11.500 3rd Qu.:28.071 3rd Qu.: 23.0000 3rd Qu.:37.392
## Max. :75.000 Max. :57.706 Max. :115.5833 Max. :68.736
## fbar recr cumf totaln
## Min. :0.0680 Min. : 140515 Min. :0.06833 Min. : 144137
## 1st Qu.:0.2270 1st Qu.: 360061 1st Qu.:0.14809 1st Qu.: 306068
## Median :0.3320 Median : 421391 Median :0.23191 Median : 539558
## Mean :0.3304 Mean : 520367 Mean :0.22981 Mean : 514973
## 3rd Qu.:0.4560 3rd Qu.: 724151 3rd Qu.:0.29803 3rd Qu.: 730351
## Max. :0.8490 Max. :1565890 Max. :0.39801 Max. :1015595
## sst sal xmonth nao
## Min. :12.77 Min. :35.40 Min. : 1.000 Min. :-4.89000
## 1st Qu.:13.63 1st Qu.:35.51 1st Qu.: 5.000 1st Qu.:-1.89000
## Median :13.86 Median :35.51 Median : 8.000 Median : 0.20000
## Mean :13.90 Mean :35.51 Mean : 7.258 Mean :-0.09236
## 3rd Qu.:14.21 3rd Qu.:35.52 3rd Qu.: 9.000 3rd Qu.: 1.63000
## Max. :14.73 Max. :35.61 Max. :12.000 Max. : 5.08000
## 'data.frame': 52582 obs. of 16 variables:
## $ X : int 0 1 2 3 4 5 6 7 8 9 ...
## $ length: num 23 22.5 25 25.5 24 22 24 23.5 22.5 22.5 ...
## $ cfin1 : num 0.0278 0.0278 0.0278 0.0278 0.0278 ...
## $ cfin2 : num 0.278 0.278 0.278 0.278 0.278 ...
## $ chel1 : num 2.47 2.47 2.47 2.47 2.47 ...
## $ chel2 : num 57.7 21.4 21.4 21.4 21.4 ...
## $ lcop1 : num 2.55 2.55 2.55 2.55 2.55 ...
## $ lcop2 : num 26.4 26.4 26.4 26.4 26.4 ...
## $ fbar : num 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 ...
## $ recr : int 482831 482831 482831 482831 482831 482831 482831 482831 482831 482831 ...
## $ cumf : num 0.306 0.306 0.306 0.306 0.306 ...
## $ totaln: num 267381 267381 267381 267381 267381 ...
## $ sst : num 14.3 14.3 14.3 14.3 14.3 ...
## $ sal : num 35.5 35.5 35.5 35.5 35.5 ...
## $ xmonth: int 7 7 7 7 7 7 7 7 7 7 ...
## $ nao : num 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 ...
Komentarz: Po pełnym czyszczeniu i uporządkowaniu
danych możemy m.in wyczytać że badany zbiór ma 52582
próbek. Poniżej prezentowane są rozkłady wartości poszczególnych
atrybutów z pominięciem X który jest numerem danej
próbki/połowu.
Poniżej zaprezentowano korelację poszczególnych atrybutów w postaci
diagramu, z wyłączeniem atrybutu X w celu wstępnego
określenia który atrybut może ale nie musi mieć
największy wpływ na długość śledzia.
Komentarz: Wstępne oględziny diagramu sugerują, że na
interesujący nas atrybut
length największy wpływ mają
atrybuty sst,nao,fbar oraz chel1. Wartości
ujemne oznaczają negatywną korelację względem danego atrybutu tj. jeśli
dany atrybut rośnie to drugi maleje, w przypadku dodatniej/pozytywnej
korelacji wzrost wartości jednego atrybutu powoduje również zwiększenie
wartości drugiego atrybutu.
Uwaga: tutaj warto zwrócić uwagę że natężenie połowów jest wyrażane w postaci ułamka pozostawionego narybku, czyli im silniejszy połów tym ten ułamek będzie mniejszy. Co może być mylące na diagramie.
Analizując powyższy diagram możemy wysnuć następujące hipotezy, które później mogą zostać poddane weryfikacji:
Może śledzie źle się czują w cieplejszej wodzie?
Może śledzie mają gorsze warunki do życia podczas specyficznych wartości oscylacji północno-atlantyckiej?
Może śledzie nie dorastają do większych rozmiarów z powodu zbyt częstych połowów?
Może ten specyficzny gatunek planktonu zawiera więcej substancji odżywczych co pozytywnie wpływa na wzrost długości śledzi?
Tworzony regresor mający za zadanie przewidzieć długość śledzi, nie powinien brać pod uwagę atrybutów które nie mają praktycznie żadnego wpływu (są bardzo bliskie wartości 0.0) na długość śledzi.
W celu “wycięcia” tych atrybutów należy wpierw ustalić tzw. parametr
odcięcia, który będzie określał kiedy dany atrybut będzie wycięty tj.
jego korelacja bezwzględna będzie mniejsza niż nasz zadany parametr
odcięcia. Parametr ten powinien być odpowiednio wysoki byśmy mogli
wyciąć najmniej skorelowane zmienne, ale tak byśmy nie wycięli zbyt dużo
atrybutów z naszego zbioru, ponieważ będzie mogło to się wiązać z
przekłamanymi wynikami. Zdecydowano się na parametr odcięcia wielkości:
0.8.
Do obliczania współczynnika korelacji wybrano metodę Pearsona.
Do utworzenia wiarygodnego modelu, zostały utworzone odpowiednio
zbiór uczący train_data oraz zbiór testowy
test_data, oba na podstawie źródłowego zbioru danych po
operacji czyszczenia o nazwie tidy_data. Na zbiór
treningowy przeznaczone zostało 75% zbioru danych. Pozostałe 25% zostało
przeznaczone na zbiór testowy.
Zbiór walidacyjny został natomiast utworzony za pomocą repeated cross-validation o pięciokrotnej liczbie powtórzeń i dwóch podziałach.
Do samego sprawdzenia poprawności modelu wykorzystany został algorytm Random Forest z miarami predykcji w postaci współczynnika determinacji (R2), pierwiastka z błędu średniokwadratowego (RMSE) oraz średniego błędu bezwzględnego (MAE).
Ostatecznie parametry do tworzenia modelu prezentują się następująco:
CutOff parameter: 0.66 - parametr
odcięcia
Correlation coefficient evaluation method:
Pearson - metoda ewaluacji dla korelacji atrybutów
Partition split: 0.75 - podział
zbioru na treningowy/testowy
Validation method:
RepeatedCrossValidation(repeats=4,number=2)
Prediction
method:Random Forest
Prediction and accuracy methods: R2, RMSE, MAE.
Atrybuty uznane za zbyt słabo skorelowane wzlędem naszego
cutOffParameter to:
## [1] "totaln" "cumf"
Komentarz: Widać że rozkłady obu zbiorów są do siebie
wystarczająco zbliżone/różne, dla naszego podziału. Nie mamy mocnego
odchylenia ani wartości odstających.
## Random Forest
##
## 39438 samples
## 12 predictor
##
## No pre-processing
## Resampling: Cross-Validated (2 fold, repeated 4 times)
## Summary of sample sizes: 19720, 19718, 19720, 19718, 19718, 19720, ...
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 2 1.170418 0.4990548 0.9256878
## 7 1.169858 0.5008801 0.9212114
## 12 1.178773 0.4943314 0.9279582
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 7.
Ostatecznie, średnia z miar prezentuje się następująco:
| x | |
|---|---|
| RMSE | 1.1582332 |
| Rsquared | 0.5086460 |
| MAE | 0.9073567 |
Komenatrz: Widać że wartość RMSE wynosi około 1, co oznacza że średnio nasz model regresora myli się co do rzeczywistych wartości śledzia o 1 jednostkę atrybutu, co w naszym przypadku jest pojedynczym centrymentrem. Jeśli spójrzymy znowu na rozkład wartości długości śledzi, to możemy dojść do stwierdzenia, że błąd wielkości około 1 centymetra jest akceptowalny.
Z drugiej strony miara R2 wynosi około 0.5 co nie jest raczej słabym wynikiem, ale leżącym w granicach normy. Wartość rzędu 0.7-0.8 byłaby tutaj o wiele bardziej satysfakcjonująca.
Komentarz: W celu określenia ważności atrybutów
użyliśmy funkcji
varImp(). Na jej podstawie jednoznacznie
możemy stwierdzić, że najważniejszym atrybutem mającym wpływ na długość
śledzi jest temperatura przy powierzchni wody (sst).
Nasze przewidywania na podstawie samego diagramu korelacji atrybutów, znajdują potwierdzenie u naszego modelu. Również natężenie połowów ma w miarę duży wpływ, tak jak przewidzieliśmy. Należy jednak wyraźnie zaznaczyć, że reszta atrybutów ma niższe wartości i istnieje wysoka szansa na zmianę kolejności ich ważności przy ponownym uruchomieniu modelu treningowego.